% Needed to run this script
% This should be put within the folder of the model
% 1. An m file for eq_c2
% 2. An m file for eq_sigma2

mode = 20001;
L    = 40000;
a    = 0.1;
freq = 5;

folder_KMC        = 'D:\2011\Projects\LPCVD\W_control\Open_loop\KMC_1\Patterned_Input\Data_2011.05.04_2\';
folder_SPDE       = 'D:\2011\Projects\LPCVD\W_control\Open_loop\SPDE_6\Data_2011.05.05\';
folder_Controller = 'D:\2011\Projects\LPCVD\W_control\Open_loop\Controller_6\Data_2011.05.05\';

W_list = 0.1:0.1:2.0;
percentage_list = [0.01,0.02,0.05,0.1];

for i = 1:length(W_list)
    W = W_list(i);
    for j = 1:length(percentage_list)
        A = W*percentage_list(j);
          
        % KMC/SPDE/Controller models are loaded from files
        file_KMC_r2 = [folder_KMC,'W_',num2str(W),'_A_',num2str(A),'\Output_r2_1.dat'];
        file_KMC_m2 = [folder_KMC,'W_',num2str(W),'_A_',num2str(A),'\Output_m2_1.dat'];
        repeat_KMC = 100;
        r2_KMC = SinkToFile(file_KMC_r2,repeat_KMC);
        m2_KMC = SinkToFile(file_KMC_m2,repeat_KMC);

        file_SPDE_r2 = [folder_SPDE,'W_',num2str(W),'_A_',num2str(A),'\Output_r2_1.dat'];
        file_SPDE_m2 = [folder_SPDE,'W_',num2str(W),'_A_',num2str(A),'\Output_r2_1.dat'];
        repeat_SPDE = 100;
        r2_SPDE = SinkToFile(file_SPDE_r2,repeat_SPDE);
        m2_SPDE = SinkToFile(file_SPDE_m2,repeat_SPDE);

        file_Controller_r2 = [folder_Controller,'W_',num2str(W),'_A_',num2str(A),'\log_r2.dat'];
        file_Controller_m2 = [folder_Controller,'W_',num2str(W),'_A_',num2str(A),'\log_m2.dat'];
        repeat_Controller = 1;
        r2_Controller = SinkToFile(file_Controller_r2,repeat_Controller);
        m2_Controller = SinkToFile(file_Controller_m2,repeat_Controller);

        % Only the fitting model is calculated
        c2 = eq_c2([W,A]);              % Need eq_c2.m
        sigma2 = eq_sigma2([W,A]);      % Need eq_sigma2.m 
        
        r2_fit = model_r2_mean_1D([c2,sigma2/c2],time,[mode,L,a,freq]);         
        m2_fit = model_m2_mean_1D([c2,sigma2/c2],time,[mode,L,a,freq]);
        
        plot(time,r2_KMC,time,time,r2_SPDE,time,r2_Controller,time,r2_fit);
        print('-djpeg',['W_',num2str(W),'_A_',num2str(A),'_r2.jpg']);
        plot(time,m2_KMC,time,time,m2_SPDE,time,m2_Controller,time,m2_fit);
        print('-djpeg',['W_',num2str(W),'_A_',num2str(A),'_m2.jpg']);        
    end
end

%%

